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ABSTRACT 

We study the secular, hierarchical three-body problem to first-order in a post-Newtonian expansion 
of General Relativity. We expand the first-order post-Newtonian Hamiltonian to leading-order in the 
ratio of the semi-major axis of the two orbits. In addition to the well-known terms that correspond to 
the GR precession of the inner and outer orbits, we find a new secular post-Newtonian interaction term 
that can affect the long-term evolution of the triple. We explore the parameter space for highly inclined 
and eccentric systems, where the Kozai-Lidov mechanism can produce large-amplitude oscillations in 
the eccentricities. The standard lore, i.e., that General Relativity effects suppress eccentricity, is only 
consistent with the parts of phase space where the General Relativity timescales are several orders 
of magnitude shorter than the secular Newtonian one. In other parts of phase space, however, post- 
Newtonian corrections combined with the three body ones, can excite eccentricities. In particular, for 
systems where the General Relativity timescale is comparable to the secular Newtonian timescales, 
the three-body interactions give rise to a resonant-like eccentricity excitation. Furthermore, for triples 
with a comparable-mass inner binary, where the eccentric Kozai-Lidov mechanism is suppressed, post- 
Newtonian corrections can further increase the eccentricity and lead to orbital flips even when the 
timescale of the former is much longer than the timescale of the secular Kozai-Lidov quadrupole 
perturbations. 



1. INTRODUCTION 

Triple stellar s ystems are believed to be very commo n 
in Nature (e.g.. iTokovhIh] 119971 : lEggleton etaD 120071 ). 
From dynamical stability arguments, these systems must 
be hierarchical triples, in which the (inner) binary is or- 
bited by a third body on a much wider orbit. Probably, 
more than 50% of the bright stars we see are at least 
(doub le) binary systems (lTokovininlll997l : lEggleton et al.l 
|2007| ). Given the selection effects against finding faint 
and distant companions, we can be reasonably confi- 
dent that the number of triple syste ms is actual l y sub - 
stantially greater than that observed. iTokovininI (|1997f ) 
showed that 40% of binary stellar systems with period 
< 10 days, in which the primary is a dwarf (0.5—1.5 Mq), 
have at least one additional companion. He also found 
that the fraction of triples and higher multiples among 
binaries with period ( 10 — 100 day) is ~ 10%. Moreover, 
iPribuUa fc Rucinskl (|2OO60 have surveyed a sample of 
contact binaries, and noted that 42±5% of 151 of them 
brighter than 10 mag. are at least triples. We can then 
conclude that many close stellar binaries with two com- 
pact objects are likely produced through triple evolution. 

Long-term stability of triple system requires hierarchi- 
cal configurations: an "inner" binary (with masses mi 
and 1712) in a nearly Keplerian orbit with semi- major axis 
(SMA) ai, and an "outer" binary in which 7713 orbits the 
center of mass of the inner binary, with SMA 02 3> ai. 



Another stability condition is that the perturber does not 
make close approaches to the inner binary orbit. In this 
stability regime, a highly inclined perturber can produce 
large-amplitude oscillations in the eccentricity and incli- 
nation o£_thes;^em^_ttie_so^alle Kozai-Lidov mecha- 
nism (lKozailll962l : lLidcwl[T96l . 

The Kozai-Lidov mechanism is an important example 
of a secular effect (i.e., coherent interaction on timescales 
long compared to the orbital period) that is common in 
hierarchical triple systems but absent from two-body dy- 
namics. This process has been proposed as a n important 
eleme n t in the evolution of triple stars (e.g.. iHarringtonl 
igp; 'Mazch fc Shaha mI IT971 iKiseleva et all | 199S 



Fabrvckv fc Trcmaina i2007t" iPerets fc Fabrvckvl l2009i 
ThompsonI 120101: iNaoz et all boilb l: iProdan fc Murravl 



2OI2I : ISharpee fc ThompsonI |2012[ ) and extrasolar plan- 



etary system s with an addition a l distant stellar com 
panio n (e.g., 'Holman et al."1997^ 'Fabrvckv & Tr cmam 
2007t iWu e t al. 2007; Takcda et al. 2008; Naoz ct aT 



2012| ) . In addition, the Kozai-Lidov mechanism has been 
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suggested to play an important role in both the growth of 
black holes (BHs) at the centers of dense stel lar clusters 
and the formation of shor t-period binary B Hs (lWenll2003l: 
iMiller fc HamiltonI 120021 ). Furthermore, llvanova et all 
(|2010( ) showed that the most important formation mech- 
anism for BH X-ray binaries in globular clusters may be 
triple-induced mass transfer in a BH-whitc dwarf binary. 

Given the hierarchical galaxy formation paradigm, and 
the strong evidence that a high abundance of the lo- 
cal galaxies host supermassive BHs (SMBHs), one ex- 
pects that major galaxy mergers should inevitably re- 
sult in the formation of SMBH binaries or multiples 
(lValtonedll996l:IHoffman fc Loebll2007tlKulkarni fc Loebl 
120121: iDotti et al.ll2012D . iBlaes et al.l (|2002D showed that 
the Kozai-Lidov mechanism plays an important role in 
the evolution of SMBH triples, where high eccentricity in- 
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duced by the outer perturber can lead to a more efficient 
merger r ate, due to gravitational wa .ve (GW) emission 
(see a lso iSetol [201^ . Also, recently lAntonini fc Peret^ 
(|2012| ) showed that secular three-body eflFects play an im- 
portant role in the evolution of binary compact objects 
near SMBH. 

GWs emitted during Kozai-Lidov-induced, highly 
eccentric orbits of compact binaries mig ht b e de- 
tectable using LIGCFI and VIR GCFI (e.g.. fV^ [20031 
but see iMandel et alT [20081 ai id lO'Lear v ct a.L' '20061). 
pulsar timing arrays (e . g.. [ Finn &: Lom mcn 201(| 
lAmaro-Seoane et al.l I2OIOI: IKocsis et alTl2012D . and fu- 
ture space-based GW observ atories, such as eLISA/NGO 
(jAmaro-Seoane et al.ll2012a| [b[j^. In fact, GWs associated 
with eccentric orbits are stronger and have a very differ- 
ent spectrum relative to their circular counterparts for 
sources at the same distance and with the same mass and 
spin. This may allow for the GW detection of eccentric 
inspirals with higher masses, larger SMAs or farther away 
from Earth relative t o their quasi-circul a r counterparts 
(lArun et al.l I2007a|[bl: [Yunes et al.l [2009t lO'Learv et all 
120091: IKocsis fc LevinI I2011D . Using GW information 
emitted by the close binary, it might be possible to con- 
strain the parameters of the third body, such as its mass 
or distance, pr ovided the GW signal-to-noise ratio is suf- 
ficien tly high (I Yunes et al.l 120111 : IGalaviz fc Briigmannl 
[20T1I) . 

The Kozai-Lidov mechanism is therefore tremendously 
important and there is still much to be understood. 
Recently, iNaoz et"al] (l2011ai rbl) showed that an eccen- 
tric outer orbit (and even a circular one with compa- 
rable mass inner binary) can behave significantly dif- 
ferently than previously assumed, the so-called "eccen- 
tric Kozai-Lidov mechanism" . Specifically, they showed 
that the inner orbits can flip from prograde to ret- 
rograde and back, and can also reach extremely high 
eccentricitie s close to unity, and the system behaves 
chaotically (jLithwick fc Naoal2011[ ). Most previous sec- 
ular three body dynamics studies that incorporated 
GR effects did so through a pseudo-potential, con- 
structed mainly to model accretion disks and 1st post- 
Newtoni an (IPN ) shi fts in the innermost stable circu- 
lar orbit (iNowak & W aTOned[T99l[ : I Artemova et al.lll996l : 
iMiller fc HamiltonI I200l ). We show here (P and Ap- 
pendix [A]) that although this pseudo-potential does cap- 
ture some IPN effects, such as the precession rate, the 
full IPN three-body Hamiltonian introduces other cor- 
rections that cannot be modeled with this potential. 

In this paper, we study the consistent inclusion of IPN 
terms in the secular dynamical evolution of hierarchical 
triple systems. We restrict attention to the IPN approx- 
imation of the three-body Hamiltonian. While it is well 
established that the eccentricity and inclination are con- 
stant in the IPN two-body problem (iDamour fc Deruelld 
|1985[) . it is not true for hierarchical triples. In addition 
to the standard GR precession of the inner and outer 
orbits, the IPN corrections lead to a new secular in- 
teraction between the inner and outer binaries that af- 
fects their long-term evolution. We find that the stan- 
dard lore, i.e., that GR effects suppress eccentricity, is 



only true when the GR timescales are several orders 
of magnitude shorter than the secular Newtonian ones. 
When the GR timescales are comparable to the secular 
Newtonian ones, we show that three-body interactions 
generally give rise to a r esonant-like eccentricity excita- 
tion (see also iFord et al.i .200 0b) . We will be using the 
term "resonance" here to describe the rapid excitations 
of the inner orbit's eccentricity, which occurs when the 
IPN timescales are comparable to the secular Newto- 
nian timescales. We demonstrate that even for systems 
with comparable inner binary masses, where the Kozai- 
Lidov mechanism is suppressed, and even when the GR 
timescales are much longer than the secular Newtonian 
ones, IPN corrections continue to excite the eccentricity. 

This paper is organized as follows. We begin with a 
definition of the parameters used to describe a hierar- 
chical triple system based on Newtonian and IPN three- 
body Hamiltonians (Sj2]). Wc then show that three-body 
evolution is modified by IPN effects (©. We discuss the 
different time-scales corresponding to the IPN effects, 
and identify the region in phase space where important 
deviations might arise due to these terms O. We then 
show that IPN terms can, in many cases, excite the ec- 
centricity of the inner orbit instead of suppressing it (Sjl)). 
We conclude with a discussion in ^ 

2. HAMILTONIAN PERTURBATION THEORY FOR 
HIERARCHICAL TRIPLE SYSTEMS 

A triple system consists of a binary (with masses mi 
and 7712) and a third body (with mass 7713) in orbit 
about the center of mass of the former. It is conve- 
nient to describe the orb its using Jacobi coordinates 
(jMurrav fc Dermottl l200dl ). Let r be the relative posi- 
tion vector from mi to m2 and R3 be the position vector 
of ma relative to the c enter of mass of th e inner binary 
(see for more details INaoz et"aL[ l2011b[ ). as shown in 

Fig.m 
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plane < 
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Fig. 1. — Coordinate system used to describe the hierarcliical 
triple system [not to scale). Here 'cm.' denotes the center of mass 
of the inner binary, containing objects of masses mi and m2. The 
separation vector r points from mi to 1712; R3 points from 'cm.' to 
ma. The angle between the vectors r and Ra is <1>. The distances 
from the bodies to a field point arc labeled by ri, r2 and ra. 
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In the PN approximation, corrections to Newtonian 
mechanics arise in powers of (f/c)", where v is the or- 
bital velocity and c is the speed of light, with n > 2 an 
integer. Here we concentrate on the IPN order correc- 
tions to Newtonian motion, which are 0(w^/c^) relatively 
smaller than the Newtonian terms. The Hamiltonian can 
then be divided into a Newtonian part ("Hn) and a IPN 
part (?^ipn): 



tot, IPN 



IPN : 



where the Newtonian part is simply 



2 TO," 2 ^ — ' r, 



1=1 



(1) 



(2) 



and the IPN part is (e.g. iSchafeii Il987t iMoord 119931: 
iLousto fc Nakajlol[2008l ): 



1 / 2 

'HiPN = — -^-^ > ^^^-i — V 



^ — ^ 



.(Pz-Pi) (n<j •Pj)(n„- -pj) 



/l^ ^ — ^ 

2? ^ 



miTrij 
mimjiTik' 



(3) 



In these equations, is the gravitational constant, ry 
(vij) is the relative position vector (magnitude) from 
mass uii to rrij, Pi (p,) is the momentum vector (mag- 
nitude) of mass rrii in an arbitrary plane (we shall later 
transform to ccnter-of-mass coordinates). In the IPN 
Hamiltonian, i,j and k' run from 1 to 3 (the three 
masses), where k' is an index while fc^ is the gravita- 
tional constant, and riij = Tij/vij. 

Many gravitational triple systems are in a hierarchical 
configuration: two objects orbit each other in a relatively 
tight inner binary while the third object is on a much 
wider orbit. If the third object is sufficiently distant, an 
analytic, perturbative approach can be used to calculate 
the evolution of the system over long timescales (rela- 
tive to the or bital period). I n the usual secular approx- 
imation (e.g.. lMarchall[T990l ). the three orbiting objects 
torque each other and exchange angular momentum, but 
not energy. Therefore, on timescales much longer than 
their orbital periods, the eccentricity and orientation can 
change, but not the SMA. 

Given this, the orbital motion of a triple system can 
be divided into two separate Kcplcrian orbits: the rel- 
ative orbit of bodies 1 and 2, and the orbit of body 3 
around the center of mass of the system. The Hamilto- 
nian for the system can then be decomposed accordingly 
into two Keplerian Hamiltonians plus a coupling term 
that describes the (weak) interaction between the two 
orbits. Let the SMAs of the inner and outer orbits be 
ai and 02, respectively. Then, the coupling term in the 
Hamiltonian can be written as a power series in th e ra- 
tio of the SMAs a = 01/02 (e.g.. lHarringtonlll968D . In 
a hierarchical system, by definition, this parameter a is 
small. 

The Newtonian part of the Hamiltonian, expanded in 



powers of a, is (e.g.. lHarringtoniri968[) . 

fc2miTO2 k^m3{mi+m2) 



2ai 



202 



(4) 



02 



00 / \ ^ / \ 

j — 2 



where Pj are Legendre polynomials, $ is the angle be- 
tween ri and r2 (see Fig. [ij and 



mi ^ — (—7712)-' ^ 
Mj = TOim2m3 ; — 



(5) 



Note that most s ecular studies follow the convention of 
iHarringtonl (|1969( ) and choose the Hamiltonian to be the 
negative of the total energy, so that "H > for bound 
systems. Here we did not follow this con vention. The 
equations of motion in iNaoz et al.l ()2011bD did use this 
convention, and thus, a reader that wishes to combine 
the two sets of equations need to introduce a minus sign 
to one of the sets. 

We adopt canonical variables, known as Dclaunay's el- 
ements, which provide a particularly convenient dynam- 
ical description of hierarchic al three-body systems (e.g. 
iValtonen fc KarttunenI [20061 ) . The coordinates are cho- 
sen to be the mean anomalies, h and I2, the arguments 
of periastron, gi and 32, and the longitudes of ascend- 
ing nodes, hi and /i2, where subscripts 1, 2 denote the 
inner and outer orbits, respectively. Their Newtonian 
conjugate momenta are 



L2- 



TOim2 



TOi + m2 
1113(7711 + TO2) 
nil +7712 + 



■\/A:2(mi + 7712)01 



(6) 



^k'^ivnx +m2 + 7773)02 , 



Gl = 



G2 = 



1 - p2 
^2 ' 



and 



Hi = Gl cos ii , H2 ^ G2 cos 72 , 



(7) 



(8) 



respectively, where ei (62) is the inner (outer) orbital ec- 
centricity. Note that Gi and G2 are also the magnitudes 
of the angular momentum vectors (Gi and G2), and Hi 
and H2 are the z-components of these vectors. The fol- 
lowing geometric relations between the momenta follow 
from the law of cosines: 



cos 7tot ^ 



Hi 



H2 = 



■-^tot 


-Gl- 


Gi 




2G1G2 




Gtot 


+ Gl- 


Gi 




2Gtot 




^2 

""tot 


+ G2- 


Gl 


2Gtot 



(9) 
(10) 

(11) 



where Gtot = Gi + G2 is the (conserved) total angu- 
lar momentum vector, and the angle between Gi and 
G2 defines the mutual inclination 7tot = 7i -|- 72- From 
Eqs. pO)) and pT|) wc find that the inclinations ii and 72 



4 



Naoz et al. 



are determined by the orbital angular momenta via 



we have defined 



cos ll : 
cos ^2 - 



'-"tot + '-^2 ~ '-"l 

2G'totG2 



(12) 
(13) 



In addition to these geometrical relations we also have 
that 



Hi + H2 = Gtot = const 



(14) 



since we are here neglecting dissipative effects such as 
GW radiation-reaction, and thus, the Hamiltonian is 
conserved. Given this parameterization, the Hamilto- 
nian or canonical equations describe the orbital motion 
via 



(15) 
(16) 
(17) 





dH 




dn 


dt 




dt 


- dL, 




dn 


dgj 


dn 


dt 


dgj ' 


dt 


dG, 


dH, 


dn 


dhj 


dn 


dt 


dhj ' 


~dt 


dH, 



where j = 1,2. 

The secular Hamiltonian (both the Newtonian and the 
IPN parts) is given by taking Eq. ([1]) expanded in pow- 
ers of a and averaging over the rapidly varying h and 
12- The averaging technique we use is know n as the 
Von Zeipel transformation (for more details, see lBrouwed 
[1959; N aoz et al. 2011b), a canonical transformation that 
eliminates the rapidly-oscillating parts of n. We apply 
this transformation twice, leading to a Hamiltonian that 
is the double average of the original Hamiltonian over 
both orbital periods. We thus refer to the resulting quan- 
tity as the "double-averaged Hamiltonian." 

The double-averaged Newtonian Hamiltonian [Eq. 
up to octupole order (i.e., up to 0{a^) beyond the leading 
order term proportional to a'^^, see INaoz et al.llMTbh 
can be written as: 



nN = n 



quad 



(18) 



where 



^quad=G2{(2 + 3e?) (3C0S2 itot - l) 

-I- 15ei sin^ itot cos(25i)} , 
15 



n 



Oct — -^f^MeiC'2{Acos(j) + lOcositot sin itot 



(1 - e^) sinpi sin 92} 



(19) 
(20) 



Here, we did not include the terms which correspond to 
the Keplerian orbital energy of the three objects which 
depend on only the SMAs, and are constant in the secular 
approximation without dissipative effects. Furthermore, 



C2 



f mi 


- "^2^ 




\mi 


+ ^2/ 





62 



(mi 4- 7712)'^ 



Li 



16 (mi + m2 -|- ms)-^ (mim2)'^ ^2^2 
5 



' 2 

B = 2 + bel-7el cos(2gi) 



and 



cos = — cos gi cos 172 — cos ztot sm gi sm 172 



(21) 
(22) 

(23) 
(24) 

(25) 



Note that the octupole coefficient in iFord et all (j2000bD , 
is simply C3 = C2(eAf/e2)15/4. Also, following our def- 
initions (see Figure [Ij mi and m2 refer to the com- 
ponent masses of the inner orbit, while 62 refers to 
the eccentricity of the outer orbit. In the test-particle 
limit (i.e., mi ^ m2) cm [Eg. (l2Tt] reduces to the oc - 
tup ole coefficient intro duced in lLithwick fc Nao3 (|2011[ ) 
and lKatz et"an (|20lH ). 



62 



(26) 



In these Hamiltonians (and in the following IPN parts), 
we have eliminated the nodes (i.e., hi and /12) by using 
the conservation of total ang ular momen t um, wh ich leads 
to /ll — /i2 = Ti"- As shown in lNaoz et al.l ()2011b[ ) this can 
be done only as long as one does not conclude th at the 
conjugate momenta are constant (e.g.. lDira3ll950[ ). The 
full equations of moti on up to th e Newto nian octupole 
order are presented in INaoz et al.l ()2011bf ). 

The averaged IPN Hamiltonian can be separated into 
different terms. First, let us use the fact that for Ke- 
plerian orbits the momentum can be related to the ra- 
dius and SMA; for the inner orbit, we can write pin = 
^in \/fc^(mi -I- m2)(2/r — 1/ai), where is the reduced 
mass of the inner orbit and a similar relation can be writ- 
ten to the outer orbit. Second, we substitute this relation 
into the three body IPN Hamiltonian, i.e., Eq. After 
transforming to the center of mass frame, the IPN cor- 
rections is expanded in powers of a up to relative 0{a^). 
This produces a similar expansion to Eq. @ for the IPN 
Hamiltonian, but due to its length we have chosen not 
to present it here. To investigate the long-term dynam- 
ics of the three-body system, we eliminate all terms with 
short-periods in the Hamiltonian, which depend on the 
rapidly changing li and I2, using a double Von Z eipel 
transformation (for more details, see lBrouweri[T959t) . In 
doing so, we must first calculate the angle between the 
vectors pin -Pout and Pout -r and pi„ -Rs, where pin (pout) 
is the momentum of the inner (outer) orbit, as defined 
in the invariable plane. 

The leading-order term in an a ^ 1 expansion is pro- 
portional to Qi^ in the double-averaged IPN Hamilto- 
nian. Keeping all terms up to 0[a^) beyond leading 
gives 



n 



IPN 



n 



IPN 



n 



IPN 



n 



IPN - 



n 



IPN 



(27) 
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where 



^IPN _ 



IPN 



IPN _ 



^IPN _ 



fc^/^in (15toi^ + 29mim2 + 15m2^ 
8ai2c2 

3fc*TOim2(TOi + TO2) 

oi^c^Vl — ei^ 
fc*mim2m3(2(mi + 7712) + Sma) 

4aia2c2(TOi + 7712 + m-3) 
k^lJ-ont (l5(7T7i + m2)2 + 29(77ii + m2)ms + Ibm^ 

8a2'^c^{mi + ni2 + ms) 
'Sk*{mi + 7772)7(73 (?77i + J77-2 + 7773) 

P 

4a2c2 (1 — e^^)^^"^ {nil + ^12) 

. aifc2777 l77l27re3 

-67773] COSltot - 



(28) 



(29) 

2\ 



(30) 

GlG2[8(777l + 7772) 

/ , \ {fniei ^fniei COS itot 

0(mi + 7772) 

+ 9ei^ (777^ + 7(7i 7772 + 7r72)cos(2gi) sin^ itot) \, (31) 



and where 



77717712 



A*out = 



7771 + 1^2 
7773(7771 + 7772) 
7771 + m2 + nis 



(32) 
(33) 



fme, = (2 - 5ef)(777f + 777^) - 3(2 - 6^)77717^2 (34) 

Here 7^^,™ includes all terms of 0{a^/'^) and 0{a^) be- 
yond , since Gi oc a^^ and G2 oc 03^^ according to 

Eqs. dll) and 0. 

Not all of the different IPN Hamiltonian terms affect 
the dynamical evolution of the triple. The H^a^-! term 
only depends on the masses and the SMAs, i.e., it does 
not depend on the canonical coordinates, and thus, it 
does not affect the canonical equations, although it does 
change _the total energy of the system. On the other 
hand, 'H^?? and Ti^^? do contribute to the dynamical 

evolution, as they clearly depend on ei and 62. A pos- 
sible, (intuitive) physical explanation for this is the fol- 
lowing. In the a2 — > 00 limit, one would expect only two 
physical effects: precession of the inner orbit and preces- 
sion of the outer orbit about the inner binary. These two 
physical effects arise because of 'H^?^ and "H^??, and 

thus, ^1™, (which satisfies > H^J^^ for large 

02/01) cannot contribute to the motion. 

The quantity Hi™ is an "interaction term," in that 
it represents the coupling between the outer and the in- 
ner orbits. Notice that this term would not be present 
if we had truncated the a ^ 1 expansion at 0{a^). No- 
tice also that the interaction term does not depend on 
the argument of periapsis of the outer orbit, 32, just 
like the quadrupole Newtonian Hamiltonian. Therefore, 
at quadrupole order, e.g. for a circular outer perturber, 
the absence of g2 in the Hamiltonian implies that the 
outer orbital angular momentun i, G2, is conserve d (t he 
so-called "happy coincidence" of lLidov"fc Ziglin.f 976f ). 

3. TRIPLE BODY EVOLUTION 



IN POST-NEWTONIAN THEORY 

The secular evolution of a three-body hierarchical 
system to Newton ian, octupole order was studied in 
iNaoz et"al] ()2011b[ ). As mentioned in ^JTJ they showed 
that the commonly assumed conservation of the z- 
component of the angular momenta of the inner and 
outer orbits {Hi and H2) is only correct in the test- 
particle approximation to quadrupole order. Newto- 
nian octupole terms further modulate the eccentricity 
and inclination oscillations. Specifically, for an eccen- 
tric and inclined outer perturber, these terms can lead 
to extremely high eccentricities and flip the inner orbit 
from prograde to retrograde. This type of behavior also 
app ears in the test-particle limit for an eccentric orbit 
/e.g.lLithwick fc Naodr2011l:lKatz et al.ll2011l:lNaoz et al.l 
1201m 

Figured] presents the secular evolution of a three-body 
hierarchical system to Newtonian, octupole order (red 
lines). We chose a system with inner binary masses 
7771 = 1 M© and 7712 = 0.001 Mq, and an outer bi- 
nary companion with mass = 10^ Mq. For this 
system, we set oi/i?f = lO'' and 02/i?f = 202, where, 
i?f = k'^{mi+m2)/c^ and i?f = k^m^/c^ are the gravita- 
tional radii of the inner and the outer orbits, respectively. 
We also set initially 62 = 0.6, ei = 0.01, gi = 92 = 0° 
and 7tot = 85°. The Newtonian quadrupole terms in- 
duce the "standard" eccentricity-inclination oscillations, 
while octupole terms modulate it. As can be seen in 
the figure, the modulation does not have a precise pe- 
riodicity and, in fact, the octupole terms introduce the 
chaotic aspects to the evolution ([Lithwick & Naoz 2011). 
When IPN corrections become significant however, the 
evolutionary orbital tracks can be significantly different 
(already at quadrupole order). 

Different IPN terms have different effects on the evo- 
lutionary orbital tracks, where the perturbations to the 
equations of motion follow from Eqs. pSHTT)) . 

"■1 

gives rise to the standard GR precession of the argument 
of periapsis of the inner orbit, while is responsi- 

ble for the precession of the argument of periapsis of the 
outer orbit. 



dgi 
dt 

dg2 
dt 



lPN(aj^^) 



lPN(a2 ^) 



3k^{mi 7772)^/2 
c2(l-e2) ' 

3fc^(777 l 4- 7772 + 7773)^/2 



5/2 



5/2 



c2(l - ei) 



(35) 
(36) 



These contributions can be recovered independently from 
the individual two-body IPN Hamiltonians of the inner 
and outer binary (see Appendix lAl Eq. IAl[) . or from an 
effective potential, or directly from the IPN metric (e.g., 
iMisner et all 119731 chapter 25 p. 668-670). Other than 
this precession, ^^-2' and , do not directly affect 

the other orbital elements^ The "H^™ term just modi- 
fies the total energy and docs not modify the long-term 
dynamical evolution at all, as long as dissipative effects 
are neglected. 
In the standard lore, if the GR precession rate of the in- 

^ Note, however, that the precession indirectly affects the evolu- 
tion of the other orbital elements through T-L^ , as shown below. 
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Fig. 2. — An example of the evolution of a triple to Newto- 
nian octupole order neglecting PN contributions (red lines) and 
including them up to the IPN interactions terms in the double- 
averaged Hamiltonian (blue lines). The system has an inner binary 
of mi = 1 M0, 7712 = 0.001 M0, and a the third object mass is: 
ma = 10"' M0. We set ai/Rf = lO"* and a2/R^ = 202. We set 
initially 62 = 0.6, ei = 0.01 gi = = 0° and itot = 85°. We 
consider, from top to bottom, the inclination of the inner orbit, 
ii, the eccentricity of the inner orbit in terms of 1 — ei, and the 
z-component of the angular momentum normalized to the total 
angular momentum. Note that we do n ot plot G2 s ince in con- 
figuration where m2 — > 0, G2 — > Ctot, IILithwick &: Naoz 2011!). 
In this case, the quadrupole Newtonian terms induce eccentricity- 
inclination oscillations, modulated by the octupole terms, while the 
IPN effect suppresses them. 



ner orbit is faster than the quadrupole secular Newtonian 
timcscales, the GR effect is presumed to suppress the ec- 
centricity growth (for an 1712 test particle to quadr upole 
Newtonian order, see iFabrvckv fc Tremaind ^OOTI) . In 
Fig- m (blue lines) we show an example where this is in- 
deed the case, even when including all IPN terms (see 
below). In this example, eccentricity (and orbital flips) 
are suppressed by the IPN corrections (the variations 
shown by the blue lines are shorter than red lines). 

The usual precession term [Eq. (j35p ] is not sufficient to 
model the system, as one must also account for the pre- 
cession of the outer orbit [Eq. and the other effects 
introduced by the IPN interaction terms (Appendix [Bjl . 
The inclusion of these terms leads to qualitatively differ- 
ent behavior because they directly drive the evolution of 
inner and outer orbital eccentricity and inclination, while 
Eq. ((35|) [Eq. ([36)) ] drives the evolution of only the argu- 
ment of periapsis of the inner (outer) orbit. Figure [3] 
shows the evolution of the eccentricity when different 
terms in the Hamiltonian are included. Wc considered a 
system with parameters mi = I Mq, m2 = I Mj (essen- 
tially a test particle), = 10^ M©, ai/Rf = 5.06 x 10^ 
corresponding to 0.005 AU. In the left column we con- 
sider initial relative inclination itot = 95° and a separa- 
tion of the outer orbit a2/R^ — 5.2 x 10"^ corresponding 
to 51.4 AU. In the right column, we consider initial rel- 
ative inclination itot = 65° and a separation of the outer 



Fig. 3. — Two examples of the time evolution of the system dom- 
inated by IPN effects, the right and left side panels differ only 
in the initial relative inclination and the outer orbital separation. 
In the top panels we consider the three-body orbit evolution due 
to Newtonian dynamics and the lowest-level IPN correction, i.e., 
Hn + W^l^. In the middle panels we add the next level of ap- 



proximation -f-'H^ 



and in the bottom panels we consider the 



approximation up to the highest level discussed here, i.e., +^1^^ ■ 
The inner binary contains an object of mass 1 Mq and an object 
of mass 1 Mj (can be considered as a test particle), while the outer 
object is a BH with mass of 10® A/q. The inner orbital separa- 
tion is ai/Rf = 5.06 X 10^, corresponding to 0.005 AU. The initial 
eccentricities are ei = 0.001 and 62 = 0.7. The arguments of 
pericenter of the inner orbit initially set to 240° and outer orbit 
initially is set to zero. In the left column, we consider an initial 
relative inclination of itot = 95°, and an outer orbital separation 
of a2/R^ = 5.2 X 10^, corresponding to 51.4 AU. In the right col- 
umn, we consider an initial relative inclination of itot = 65° , and 
an outer orbital separation of a2/R^ = 4.8 X 10^, corresponding to 
47.35 AU. Observe that as one includes more IPN effects, qualita- 
tively different behavior emerges. 



orbit a2/Rl = 4.8 x 10^ corresponding to 47.35 AU. For 
the calculation in the two columns, the initial eccentrici- 
ties are ei = 0.001 and 62 = 0.7, and the initial argument 
of pericenter of the inner and outer orbits is set to 240° 
and zero, respectively. This system configuration is such 
that the lPN(aj~^) timcscales for circular orbits (^ 59 yr, 
for the left column example, sec 21 for niore details) are 
shorter than the Newtonian quadrupole ones (^ 145 yr 
for the left column example). Nevertheless, the secu- 
lar eccentricity oscillations are still present. The bottom 
panel shows that the interaction term V.}^^ introduces 
a significant new periodic modulation to the eccentricity 
evolution. We discuss in more details the conditions in 
parameter space that lead to this behavior in 21 ^^'^ El 
(note that this system represents the resonance peak of 
the 95° and 65° cases of Fig. E] below) . 

One might worry that a IPN treatment might not be 
sufficient to model certain regions of phase space, as we 
have neglected 2PN and higher PN order terms in the 
evolution. Such terms become important when the PN 
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perturbation parameter, vjc with v any system velocity, 
i.e. the pericenter velocities {vpi/cY = )/[ai(l — 

ei)] and similarly for the outer orbit, are not sufBciently 
small. In most of our examples we ensured that our 
systems are well within the PN regime, however, for very 
eccentric systems, 2PN and higher PN order terms may 
be important. In fact, one might naively expect the IPN 
corrections accounted for here to be negligible if u/c ^ 1. 
This is not so, because although the IPN terms are small 
at any point in time, their effects may accumulate and 
become significant over long timescales in the three-body 
problem. 



4. TIMESCALES 

In order to explore the regions of phase space where 
the IPN effects may be expected to become significant, 
we compare the various timescales associated with the 
individual terms in the Hamiltonian. 

The timescalc associated with the Newtonian 
quadrupole term can be estimated from the canonical 
equations. More precisely, iq^^d ^ 27r/C2, where C2 is 
given in Eq. p2|): 



which gives 



^quad 



3/2 



(37) 



5/2 2/'i 2\ 

t -2 ^ Ztt 



t 



"1 3fc'^(TOi + 7712)'^/^ ' 

IPN o '^2 '^^(1 ~ ^1) 

- ' ZTT- 



t 



3fc3(mi + ?7i2 + ms)^/^ 
iPN_ 16 aie^(l-ei)3/2(mi+m2)3/^ 

9 ei-\/l — e1k^{ml + mim2 + m|)m3 



(39) 
(40) 
(41) 



All of these timescales depend on the secularly varying, 
inner orbital eccentricity, which implies that we need 
to explore different eccentricity values in phase space. 
Equations ([55HiT|) show that the Newtonian and IPN 
timescales have a simple dependence on the inner and 
outer orbital eccentricity, on the mass ratio qm ~ my, /mi, 
and on the SMAs. 

If any of the above timescales is much smaller than 
all others, then the evolution of the triple is dominated 
by the corresponding term in the Hamiltonian. Next, 
we examine the three-body parameters where any two 
timescales are equal, which defines the region where the 
corresponding two terms are equally important. The cor- 
responding regions are shown in Figure 2] in the m2 — > 
limit. 

Equating t^?? to t^^^^ (Eqs.inZlandUn]) gives a relation 

between the SMAs which normalized to the gravitational 
radius of the inner and outer binaries, R\ and as 
defined above, can be written as: 



The timescales associated with the Newtonian oc- 
tupole terms are more difficult to estimate due to their 
chaotic effect on the orbits. For example, in Fig. [5] 
the first modulation p eriod is shorter than the second 
(see. lNaoz et al.ll2011bl for more example s). However, as 
demonstrated in iLithwick fc Nao3 ([20Tl[ ). although the 
system is chaotic when the octupole terms are included, 
there are two general features in the evolution: one as- 
sociated with an octupole timescale, defined below, and 
a shorter one (see for example their Fig. 7 of surfaces of 
section, which shows the two different evolutionary be- 
haviors). We define a timescale for the regular part of 
the Newtonian octupole evolution through the rough es- 
timate t^^^ ^ (4/15)iquaci/eM for a given inner and outer 
eccentricity. 



2-K 



4 a|(l - ei)5/2yT3^(mi + ^3)'/' 



15 



5/2 



e2k\mi 



(38) 



TO2 ms 



Note that when the inner binary is very eccentric, these 
timescales can change by orders of magnitude. At oc- 
tupole order, the eccentricity of the outer orbit can also 
oscillate, although usually these oscillations are small in 
magnitude. For the remainder of this section, we employ 
a test-particle approximation for one of the components 
of the inner binary, i.e. m2 — 0, for which 62 = const. 
(jLithwick fc NaozlBoTll) . We will see that this is suf- 
ficient to understand the regions of phase space where 
IPN terms become important. 
The IPN timescale can be estimated as in Eq. (|37)). 



02 



quad=lPN(aj^^) 



1/3 



ai 



4/3 



1 (l-ej)i/3 
2/3 



(42) 

This relation is shown by the blue hatched area bounded 
by solid blue lines in Figure S] for < ei < 0.999. A 
resonant-like IPN excitation of eccentricity is possible 
in this region, as we will show in the next section. For 
much larger 02 or smaller ai, t^Jl^ ^ ^quad' ^'^^ ^^^^ 
the Kozai-Lidov eccentricity excitations are suppressed 
by the IPN effects. 

^quad gives 



Next, equating t^^^ to 



Q2 



quad=lPN(a2 



1 



(1 



1 



(43) 



This is shown by a blue short-long dashed line in Fig. |4] 
on the top and bottom panels for 62 = 0.5 and 0.7, re- 
spectively. For a2/i?f much larger than this value, the 
Kozai-Lidov oscillations are suppressed and the IPN ef- 
fects dominate. 

Let us next compare tf^^ and t^^^ by setting them 
equal to each other: 



oct=lPN(int) 



32 
1357r 



eie2 



g,„(l-el)(l-e?) 



(44) 

shown by long-dashed brown lines in Fig. 21 The IPN 
effects are equally important as the Newtonian ones in 
the gray shaded area in Fig.|3]for < ei < 0.999. This is 
the regime in which the IPN interaction term introduces 
qualitatively different behavior in the orbital motion (i.e.. 
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Fig. 4. — Regions of binary parameters where the various secular Newtonian and IPN effects are expected to become significant. We 
show the SMAs where the timescales, corresponding to individual terms in the Hamiltonian, are equal to each other for different ai/Rf 
and 02/^3. We consider four mass ratios between the outer object and the inner binary qm = ma/mi = 0.01, 1, 10* and 10^ (see labels 
in each panel) and two different choices for outer orbital eccentricity (e2 = 0.5, top panels and £2 = 0.7 lower panels), setting m2 — > in 
all panels. We consider t\^^ = t^™ (Eq. gg, solid blue lines), t\^^ = t'^™ (Eq. gSJ, thick long-short dashed blue lines), t^^ = ^IPN 

(Eq. 1144 l l. long dashed brown lines) and t^^.^ = (Eq. 1 145 I I. short dashed purple lines). For the latter, we show ei = (top purple 

^2 

line) in all panels and in the top left panel we also show the eccentric case ei = 0.999 (bottom purple line). The gray and blue shaded 
areas cover the range < ei < 0.999 between the brown shaded lines and the blue solid lines respectively In the top left panel and bottom 
right, we also show a green band for which the timescale to shrink the inner orbit by a factor two, due to GW emission, is equal to the 
quadrupole timescale, which covers the range < ei < 0.999. This line is generated by specifying mi = 100 Mq and m2 = = 1 Mq 
(top left panel) and m2 = mi = 1 M0 ( bott om right panel). We also show two stability criteria, e = 0.1, dotted magenta lines and the 
IMardling fc AarsethI 11200 II ) criterion [Eq. 1 146 I I , in the top right panel] . The green dot in the right top panel represents (up to a factor 3) the 
location in this phase space of the system considered in Fig. [7] while the black dot in the bottom right panel represent the location of the 
example considered in Fig. [3] The rectangle in the top left panel shows the parameter space considered in Figure |6] The black rectangles 
in the top and bottom right panels roughly represent the region in parameter space where planets in stellar binaries (top right) and stars 
in the galactic nuclei (bottom right) live. 



modulation) as shown in Fig. O Outside the gray region 
in Fig. m the interaction term is neghgible. Note that 
comparing tf^^ with tq^^^j leads to a vertical line in the 
phase diagram of Fig. This is because both timescales 
have the same dependence on the outer SMA (~ 02), 
resulting in ai/Ri < 10 (not shown). In ^JSJ we explore 
the parameter space that also leads to excitations in the 
eccentricity (as shown in Fig. [S]). 

Comparing t^^^ to the octupole timescale gives usually 

a longer timescale than the quadrupole (short-dashed 



purple lines in Fig. 2]): 



02 



oct=lPN(a2 ) 



1/3 



5/3 



(45) 



2/3 



W (l-el)(l 



.2)1/3 



This relation also provides a range of parameters for dif- 
ferent values of ei and 62, but to avoid cluttering, we plot 
this timescale only for a circular inner binary (ei = 0). 
The boundary shifts to larger 02 in the eccentric case. In 
top-left panel of Fig. U (the qm = 0.01 case), we show the 
range of Eq. ((45|) for 62 G (0,0.999), since, in this case, 
this ratio is smaller than Eq. (|44]). 
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Figure 2] also shows the regime of validity of the hierar- 
chical triple approximation, where we choose e > 0.1 as 
a rule of thumb for stability (dotted magenta line). For 
the Qm. = 1 case, this rule of th umb seems to agree with 
the lMardling fc AarsethI (j2001l ) stability criterion, which 
defines a stable three-body system as one that obeys 

->^-^(' + 'l-y (l-e.)^/4 ^'l80^j ' 

where in the top right panel we considered this crite- 
rion with prograde itot = 0°. Note that retrograde mo- 
tions are usually mo re stable (e.g. Ilnnaneni [1971 [l980l: 
IMorais fc GiuppondlMl) . 

Another consistency requirement for the IPN Kozai- 
Lidov effects to operate is that gravitational radiation 
reaction does not change the SMA significantly over this 
timcscale. We define tcwi, the GW inspiral timescalc 
of the inner binary Jot the SM A to chang e by factor of 
two using iPcters ( 196i ) fsee also lArun et al . 20091. Note 
that icwi — >■ oo when m2 — )- 0, but the GW inspiral may 
become very important in the comparable-mass limit and 
shut off the Kozai-Lidov effect. This is the case above 
the green bands in the top left and bottom right panels, 
which show the range of SMAs where iqu^d = *GWi for 
7712 = iTT-s = 1 Mq (top left panel) and 7712 = mi = 1 Mq 
(bottom right panel) for < d < 0.999. Note that 
there is a region below the green band where the IPN 
approximation is insufficient, and where 2PN and higher 
PN order corrections need to be included; we leave this 
to future work. 

We conclude that IPN effects may be important for 
a vast range of parameters as shown in Fig. 2] Note 
however, that physical timescales are not shown there; 
this figure is independent of an overall dimensional scale 
(e.g. total mass or the scale of the SMAs). The physical 
timcscale may be smaller or larger than the Hubble time 
depending on the actual system parameters considered. 

5. EXCITATION OF THE INNER ORBITAL ECCENTRICITY 

As stated above the IPN effects may suppress the 
Newtonian Kozai-Lidov oscillations if the correspond- 
ing IPN timescale is much smaller than the Newtonian 
quadrupole one. However, we identify two regimes where 
the combined secular Newtonian and IPN effects pro- 
duce interesting different behavior in three-body sys- 
tems: (i) if the quadrupole order terms are comparable 
to or somewhat smaller than the IPN ones, or (ii) for 
comparable-mass inner binaries where the IPN effects 
are subdominant relative to the Newtonian quadrupole 
terms but they are comparable to the Newtonian oc- 
tupole timescale. We discuss these two regimes in detail 
below. 

5.1. Eccentricity Peak in the Restricted Three Body 
Problem for a Massive Perturber 

The standard lore says that GR effects (or IPN ef- 
fects in our case) suppress the eccentricit y growth of 
the inner orbit in a three body sy stem fe.g.. IBlaes et al ' 



120021: iFabrvckv fc Tremainell2b07l) . However . iFord et al 



(l2000bO, studying the triple system PSR B1620-26 



and extend the investigations of that study for a broader 
range of three-body systems and examine (i) whether 
the new IPN terms derived in fj2] give rise to different 
behavior; and (ii) whether Newtonian octupole terms 
can significantly change the evolution of the three-body 
system in the presence of the IPN terms by producing 
flips of the inner o rbit and eccentricities close to unity 
()Naoz et al.ll2011ai rbl). We begin by exploring systems in 
which 7Ti3 ^ TTii {qm ^ 1), and investigate the opposite 
limit in the next subsection. 

Let us systematically examine the parameters where 
the competition between the secular Newtonian Hamil- 
tonian and IPN corrections lead to the excitation of the 
inner orbital eccentricity as opposed to a suppression. 
We do this by preforming a large number of simulations 
starting from a nearly zero eccentricity for the inner bi- 
nary and varying the following dimensionlcss parameter: 



7^ = 



t^ 

quad 



3 



1)3/2 



(47) 



showed that the Newtonian octupole and the leading or- 
der IPN corrections [i.e., Eq. ((35)) ] can lead to the exci- 
tation of the eccentricity of the inner orbit. We repeat 



where we have substituted Eqs. ([ST]) and (|39l) with ei = 
0. The quantity in Eq. (|T7| compares the timescales 
of the leading order IPN and secular Newtonian effects. 
We find that a resonant-like eccentricity excitation can 
take place around TZ = 1 for zero initial eccentricity as 
shown below. Note, however, that another important in- 
gredient to this resonant-like behavior is the Newtonian 
octupole term, which is most obvious in simulations with 
a low initial mutual inclination, for which the quadrupole 
approximation is subdominant (see also ij5.2l below and 
iFord et ani2000b() . In other words, neglecting the contri- 
bution of the octupole-level of approximation, one would 
miss entirely the resonant behavior. 

Figure [5] shows the maximum eccentricity achieved 
during the course of the evolution of the system studied, 
after 1000 quadrupole cyclcfQas a function of TZ. We also 
examined the effects of_various IPN terms by repeating 

the calculations using Hn + H'^^^ (blue triangles), then 

°i _ 

adding H^?^ (green crosses), and finally including Hf^^ 

(red circles). The fiducial example chosen (black lines 
in all panels) has mi = 1 Mq, ?ti2 = 0.001 Mq, and 
7713 = 10^ Mq (i.e. Qm = 10^) with initial conditions 
ei = 0.001, 62 = 0.7, itot = 65°, gi = 240° and 52 = 0°. 
To generate Fig. [5l we used different runs and changed 
both ai and 02. For a choice of mass ratio and initial 
mutual inclination the parameter TZ collapses the differ- 
ent runs with different SMAs into one curve. Note that 
changing the initial value of the argument of periapsis 
does not change the location of the resonance, but it does 
change its width, as shown in the inset of Fig.[5l where we 
compare the fiducial example with an initial gi = 240° 
to an initial gi = 0. The different panels show that the 
amplitude and width of the resonance are different for 
different qm and £2, as they depend on the Newtonian 
octupole terms. The top left panel shows that different 

* Each point in this figure corresponds to a separate, high- 
resolution three-body evolution with IPN effects, each of which 
takes approximately 3 days to complete per computer core. We 
also preformed convergence tests using longer integration times at 
high resolutions, and found that for this type of systems at least 
1000 quadrupole cycles are needed for convergence, over the pa- 
rameter range considered in the Figure. 
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Fig. 5. — The maximum eccentricity in a triple system in the test-particle approximation for different mass ratios (gm = 10*', left hand 
panels, and = 10*, right hand panels) as a function of the relative timescales of the IPN and secular Newtonian quadrupole effects 
(7?., see Eg. I47I I. We consider the IPN evolution including terms only up to 0(aj~^), 0{a^^) and the interaction term (blue triangles, 
green crosses, and red squares, respectively). We show three different initial outer orbital eccentricities: 62 = 0.7 (top panels) 62 = 0.6 
and 0.5, bottom right and left hand respectively. In all these examples, we set mi = 1 Mq, m2 = 0.001 Mq, ms = qmTni and we vary 
both the inner and outer SMAs to match the different values of TZ, see Eq. I I47I I. The system is initialized with 6i = 0.001, g2 = 0° and 
gi = 240° (inset shows the results for initializing gi = 0° and 240°). In all panels we set the initial mutual inclination to 65°. In the 
top left panel, we also set the initial inclination to 95° (left set of lines) and 50° (right set of lines) . We compare the fiducial example 
ilm = 10^,itot = 65°, 62 = 0.7), black lines in each panel, to systems with a different mass ratio (g^ = 10''), right hand panels), and 
different eccentricities, see labels top right panel. Observe the emergence of resonant-like eccentricity excitations. 



initial inclinations change both the location of the reso- 
nance and its width (as a function of 7?.), i.e., low initial 
inclination leads to broad resonances. The time evolu- 
tion of the systems with parameters on the peak of the 
resonance of the itot = 95° and 65° cases are shown in 
Fig. El 

To the left of the resonance peak, the maximum ec- 
centricity converges asymptotically to the initial inner 
eccentricity for TZ <^ 1 (i.e. in this case the eccentric- 
ity is not excited). This was chosen to be ei = 10~^ in 
Fig. [5] Decreasing the initial eccentricity (not shown) 
changes the figure at 7?, ^ 1, but it does not change the 
amplitude and location of the resonance and the behav- 
ior at larger TZ. To the right of the resonance {TZ > 1), 
the maximum eccentricity starts to increase when the ec- 
centric Kozai-Lidov mechanism begins to dominate. For 
9m 2> 1 and 7?, 3> 1, we find asymptotically ei — 1 (see 
^5.21 for small q™). 

A binary that slowly shrinks due to GW emission, 
slowly changes TZ and may sweep across the resonant re- 
gions of eccentricity excitation shown in Fig. [5l Thus, 
resonant IPN eccentricity increase may take place in 
many inspiraling systems some time during their evolu- 
tion. Whether this occurs or not depends on the masses 
and SMAs; a detailed analysis is left to future work. 



5.2. The Case of PSR B 1620- 26 

Although the above discussion of the resonant behav- 
ior assumed a test particle (m2 ~ 0), general mass-ratio 
triples also exhibit a similar effect (as was first consid- 
ered bv lFord et al.l[2000b[ ). provided t^™ is shorter then 
"i 

^quad [^6^ Fq. (p9| ]. To find TZ in the general case, we 
can simply set 

qra = rnzl {nil + "12) (48) 

in Eq. (gT]). iFord et all (|2000b( ) observed a resonant- 
like eccentricity increase while studying the triple sys- 
tem PSR B1620— 26, which is located near the core of 
the globular cluster M4. They showed that this resonant 
behavior may explain the unusually large eccentricity of 
the inner binary, which contains a millisecond radio pul- 
sar of m^ = lAMc?, and a companion of 1712 = 0.3Mq 
(IMcKenna fc Lvnel[l988l) . 

For c ompleteness , we rep eat and extend the calcula- 
tion of iFord et all ()2000bl ) by fixing the inner orbit's 
SMA, changing the outer orbit's SMA, and including all 
IPN terms. We choose two different values for the in- 
ner binary to e xplore the sensitivit y to these parameters. 
First, following iFbrd et all (|2000a[) . we consider ai/Ri = 
5.6 X 10^ (i.e. ai = 0.77 AU). Additionally, we consider 
ai/Ri = 3.6 X 10^ (i.e. oi = 5 AU). We adopt param- 
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eters for the outer perturber from iFord et all (|2000aD : 
7713 = 0.01 Mq and 62 = 0.45. We initialize the system 
with gi — 92 = 0°, ei = 10"'^, and we also choose two 
different initial inclinations, itot = 65° and itot = 0°. A 
mutual inclination of i^ot = 0° highlights that the per- 
turbations of the outer orbit affect the inner orbit due to 
the Newtonian oetupole term, even far from the nomi- 
nal Kozai-Lidov regime. In this configuration, the "H^™ 

term is the most significant IPN effect, as shown in Fig- 
ure m (see the rectangle in the top left panel) . 

The left-panel of Fig. |6] shows the inner orbit's maxi- 
mum eccentricity as a function of the outer orbit's SMA 
(or equivalently the TZ value for oi = 0.77 AU). This 
figure co nfirms the resonant -like increase in eccentricity 
found in lFord et all (|2000bf ) (their Figure 14). This fig- 
ure also shows that the resonant-like eccentricity increase 
is present for a large range of 02 values, even for sys- 
tems where the excitation of the eccentricity due to the 
Newtonian oetupole term is somewhat suppre ssed, due to 
compa rable masses for the inner orbit (e.g.. iNaoz et all 
l2011bD . Ahhough chan ges in the outer orbit's eccentric- 
ity (see the right panels of Fig. |6]) do not change the loca- 
tion of the resonant peaks, their amplitude does change 
for large inclinations. The ai — 5 AU case has irregular 
behavior and results in a higher inner and outer orbital 
eccentricity. Note that "H^^f term is the dominant one 

here and the other IPN terms are negligible in this con- 
figuration, as can be see from the black rectangle in the 
top left panel in Figure 21 

The right panels of Fig. [6] show the time evolution 
of the inner and outer orbital eccentricity for the ai = 
0.77 AU case with itot = 0° and a2/R3 = 4.9 x 10" 
(a2 = 48 AU, bottom panel) and itot = 65° and 02/ R3 = 
2.2 X 10^^ (a2 = 22 AU, top panel). Unlike the systems 
considered in Fig. [SJ the outer orbit's eccentricity oscil- 
lates slightly (see right hand panels in Figure [6|) . Fur- 
thermore, the eccentricity peak is even larger than the 
eccentricity reached for TZ > 1, which shows that the 
IPN terms further increase the inner orbit's eccentricity 
above the excitation induced by the Newtonian eccen- 
tric Kozai-Lidov mechanism. This was not the case in 
Figure [5] primarily because (?m ^ 1 there. 

Although we integrated the system for up to 1000 
Kozai-Lidov cycles, the eccentric Kozai-Lidov process did 
not seem to induce chaotic behavior in this configura- 
tion (since 7713 < nii + 1712). Numerical convergence was 
reached already after a few hundred Kozai-Lidov cycles 
of evolution (in thin black lines we show the results of 
integrating the system up to 1 Kozai-Lidov cycle). The 
bottom right panel shows that integrating over only one 
Kozai-Lidov cycle, for the itot = 0° case, misses the long 
timescale oscillation. This explains why ou r eccentric- 
ity pea ks are slightly higher than those of iFord et ahl 
(|2000bD . This also explains t he somewhat larger ei^max 
values we found compared to iFord et ahl ()2000bD in the 
regime where the IPN effects are subdominant and the 
eccentricity of the inner orbit increases due to the eccen- 
tric Kozai-Lidov mechanism. Note however, that the life- 
time of this system in the core of the globular M4 is about 
one Ko zai-Lidov cycle, which explains why iFord et al.l 
()2000b[) did not bother to evolve over many Kozai-Lidov 
cycles. 
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Fig. 6. — The maximum eccentricity as a function of the outer 
orbit's SMA (left panel) in a triple for the case of PSR B1620-26. 
The inner binary is a miUisecond pulsar of mass 1.4 Mq with a 
companion of m2 = O.SMq, and the outer body has mass 7713 = 
O.O1M0. The inner orbit has ai = 5AU and 0.77 AU, in two 
different sets of simulations (see labels) . The initial eccentricities 
are ei = 10~* and ea = 0.45 and the initial relative inclination 
*tot = 65° and itot = 0°. The argument of pericenter of the inner 
and outer orbits are initially set to zero. The top axi s sh ow the 
value of Ti for the case of ai = 0.77 AU, see equation I I47II . Note 
that the itot = 0° resonance happens when 7?. ~ 1. The thin 
black line shows the results of integrating the system over a single 
Kozai-Lidov cycle. Although we include all IPN terms described 
in previous sections, the curves corresponding to the lower order 
ones exactly overlap the curves including the interaction term (i.e., 
the leading IPN term is the most dominant in the evolution of the 
system). In the right panels, we show the time evolution of the 
eccentricity of the inner (blue line) and outer binary (green line). 
Here we set ai = 0.77 AU itot = 0°, 02 = 48 AU in the bottom and 
7tot = 65°, a2 = 22 AU in the top panel, respectively. In this case, 
the resonant eccentricity excitation due to the IPN terms reaches 
higher values then the one achieved by the eccentric Kozai-Lidov 
mechanism for small (large) a2 (7?.). 



5.3. Orbital Flips and Eccentricity Excitation for 
Comparable-mass Inner Binary 

If the inner binary consists of comparable mass- 
objects, the Newtoni an oetupole term is sup pressed (see 
Eq. m]). Recently, iSharpee fc ThompsonI (|2012f ) con- 
sidered the evolution of triple systems with comparable 
masses and showed that the eccentric Kozai-Lidov evolu- 
tion can be triggered if one of the stars in the inner binary 
loses mass. We show here that the eccentric Kozai-Lidov 
evolution can also be triggered without mass loss, but 
accounting for IPN effects, as shown in Fig. [71 For this 
figure, we set 7711 = 10 M©, 7773 = 30 Mq, ai/i?f = 1 x 10^ 
(corresponding to 10 AU), and a2/i?f = 1.7 x 10^ (cor- 
responding to 502 AU). The initial eccentricities were 
ei = 0.001 and 62 = 0.7 and the initial relative incli- 
nation 7tot = 94°. The argument of pericenter of the 
inner and outer orbits was initially set to 240° and zero, 
respectively. The two panels in this figure differ in the 
choice of 7712, i.e. in the left panel 7772 = 0.001 and in the 
right panel 7712 = 8Mq. While a test-particle evolution 
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Fig. 7. — Eccentricity excitation and orbital flips for the New- 
tonian octupole and the IPN approximations (up to the interac- 
tion term, see text for details) as a function of time. We com- 
pare the test-particle case (m2 — > left panels) to a comparable 
mass case (m2 = 8 right panels) in the inner binary with 

mi = 10 Mq, always with an outer object of mass = 30A'/q. 
The separation of the inner orbit is ai/R^ = 1 X 10* (correspond- 
ing to 10 AU and orbital period Pi ~ 10 yr), and the outer orbit's 
separation is a2/R^ = 1.69 X lO'^* (corresponding to 502 AU and 
P2 = 2.8 X lO'^yr). The initial eccentricities are ei = 0.001 and 
62 = 0.7 and the initial relative inclination itot = 94°. The argu- 
ments of pericenter of the inner and outer orbits are initially set to 
240° and zero respectively. For these examples TZ ^ I. We show 
with red lines evolutions without IPN corrections (curves including 
only the lower order IPN terms simply overlap this). We also show 
the minimum eccentricity corresponding to the detectable LIGO 
frequency range (horizontal lines in the bottom panels). The IPN 
corrections help to further increase the eccentricity and lead to 
orbital flips for the inner binary for comparable masses. 



is relatively insensitive to the IPN terms in this case, 
comparable mass systems present qualitatively different 
behavior. In particular, while the Newtonian eccentric 
Kozai-Lidov effect is suppressed for comparable masses, 
IPN effects can trigger it. A possible reason for the qual- 
itative difference is that changing m2 from zero to 8 Mq 
resulted in a configuration for which the leading order 
IPN timescale is closer to (but still slightly longer than) 
the octupole time scale. In the test particle case, the oc- 
tupole timescale is two orders of magnitude shorter than 
the shortest IPN timescale. 

The excitation of the eccentricity depends on the im- 
portance of the Newtonian octupole term. Without 
IPN effects, comparable mass triples result in ej\/ — > 
[Eq. (HI])], which suppresses the eccentricity excitations 
and the flipping of the inner orbits. This can be seen 
in Fig. |8] (black and gray lines) for two examples, where 
we consider Am = \ini ~ m2\ = 2 M©, (black lines) 
and Am = 0.2 Mq (gray lines). For both of these, 
7713 = 3m 1 and the two systems were initialized with 
62 = 0.7, ei — 0.001 and itot = 94°. However, we find 
that although IPN effects are small during a single orbit, 
they can be significant on much longer timescales, and 



can lead to significant eccentricity growth and orbital- 
flipf0. The top panel shows that for tM > 0.002 a flip 
is triggered due to the IPN terms. In other words, for 
Cm > 0.002 the colored curves deviate (including IPN 
effects) from the black and grey ones (Newtonian effects 
only). For larger em, where the eccentric Kozai-Lidov 
mechanism is triggered, the maximum eccentricity can 
be very close to unity, and thus due to the chaotic nature 
of the system, the maximum eccentricity shown should 
be considered as a lower limit. For cm ^ 0.005 we could 
not reach convergence after 1000 quadrupole cycles, since 
the eccentricity is very close to unity (see Teyssandier et 
al; in preparation). 

In the above examples, a requirement for eccentricity 
excitation is that the IPN timescale, t^^^ , be shorter 

than or comparable to the octupole timescale t^ct^ 
Eq. (I25|)], i.e., t^^^^ < tl^^ ~ t^^^^. A possible explana- 
tion for the excitation of the eccentricity in these cases 
is the following. Neglecting IPN effects, comparable- 
masses in the inner binary suppress the Newtonian 
octupole effects and the outer potential is effectively 
quadrupolar. However, GR precession of the inner orbit 
breaks this symmetry. As long as GR precession occurs 
on a timescale comparable to (or slightly smaller than) 
the octupole one, the eccentric Kozai-Lidov mechanism 
is then triggered. 

Let us now discuss the implications of these finding for 
direct GW detections using Earth-based instruments like 
LIGO and VIRGO. The characteristic frequency of the 
GW signal is fp = Vp/vp, where Vp and are th e orbital 
velocity and radius at pericenter (|Peterj[l964[ ). Thus, 
fp = 27r(l + e)i/2(l - e)-3/2p-i, where P is the orbital 
period. We assume that the GW signal is in the de- 
tectable frequency band ii fp > 5 Hz. For a comparable- 
mass inner binary, the Newtonian eccentric Kozai-Lidov 
mechanism is suppressecP^ and the eccentricity remains 
smaller than in the test particle case. The GW frequency 
emitted by a circular binary with an orbital period larger 
than a second is too small for a LIGO detection. How- 
ever, the IPN eccentricity excitations discussed in this 
paper lead to a much larger fp and might lead to GWs 
in the LIGO band. In particular. Figs. [7] and [8] show 
examples where 1 — ei can be as small as 10~^ to 10~^ 
for a comparable-mass inner binary. These sources enter 
the LIGO GW frequency band if their orbital period is 
less than 1 to 60 years. However, note that this estimate 
is oversimplified because it neglects the backreaction of 
GW emission on the evolution. The latter strongly re- 
duces the SMA of the binaries during close approaches, 
and may lead to an eccentric inspiral and merger within 
a Kozai-Lidov period. The GW inspiral may deliver the 
binaries to the LIGO frequency band even if the signal is 
outside the LIGO band during the Kozai-Lidov oscilla- 
tions. If the event rate of these sources is sufficiently 
large within the LIGO detection range, these so urces 
could constitute a distinct po pulation for LIGO (jWenI 
[2003tlAntonini fc PeretsllMa) . 

^ The choice of argument of pcriapsis does not change the out- 
come; the first example was initialized with gi = 240° while the 
second had gi = 0° , both had 92 = 0° . 

Although suppressed, note that the eccentricity can still reach 
0.999 in this case. 



Resonant IPN Eccentricity Excitation in Hierarchical Triples 



13 




Fig. 8. — Excitation of the eccentricity and inclination due to IPN effects for Kozai-Lidov timcscales much shorter than the lowest IPN 
timescales (as parametrized through ejv/i see Eg. 12111 . We examine the maximum eccentricity reached, and plot 1 — ei ^iax (bottom panel), 
and the maximum and minimum inclination reached in the runs (top panel). We show two examples. The first system is initialized with 
nil = 10 M0, m2 = 8 Mq ms = 30 M@ and 62 = 0.7. We initialize the system with ei = 0.001, gi = 240°, 32 = and itot = 94°. We 
vary both the inner and outer SMAs to match the different values of ej\/ depicted in the figure. We consider the Newtonian approximation 
(filled black squares) and the three IPN level of approximations, 0{a^^), ©(a^'^) and the interaction term (filled blue triangles, green 
cross and filled red squares, respectively). In the second system, we set mi = 1 Mq, m2 = 1.2 Mq 1713 = 3 Mq and 62 = 0.7. We initialize 
the system with ei = 0.001, gi = 0°, 92 = and itot = 94°. We vary both the inner and outer SMAs to match the different values of e]\f 
depicted in the figure. For this example, we consider the Newtonian approximation (empty gray squares) and up to the interaction level 
of the IPN approximation (empty magenta circles). We also show the detectable LIGO frequency limit for the first example, where we set 
a2 = 501 AU and varied ai between 5 AU to 15 AU (solid brown line). 



6. DISCUSSION 

The Kozai-Lidov mechanism (|Kozailll962l lLidovlll962L 
see below), has been shown to play an important role for 
highly inclined hierarchical triples, from planetary sys- 
tems to stellar size and/or massive compact objects (e.g., 
iNaoz et"ani2011bl and references therein). For an eccen- 
tric outer perturber, the eccentricity of the inner orbit 
can reach values extremely close to unity, a nd the incli- 
nation can flip from prograde to retrograde (jNaoz et al.l 
l2011allH ). The quadrupole Kozai-Lidov oscillations be- 
tween the eccentricity and inclination still persist at oc- 
tupole order, but they are further modulated on long 
timcscales. 

We have here studied how the Kozai-Lidov mechanism 
is affected by IPN corrections to the three-body Hamilto- 
nian, focusing on secular and hierarchical three body sys- 
tems. We expanded the IPN Hamiltonian in the ratio of 
SMAs (a) to third order beyond leading, i.e. the leading- 
order terms in the IPN Hamiltonian perturbation scale 
here as a^^ and we carried out an expansion up to rela- 
tive 0{a^). We also averaged over the orbital timcscale 
of the inner and outer binary to investigate the long- 
term secular evolution of the system (^J^I- We examined 
the effects of the different IPN terms in this expansion: 
IPN precession of the inner orbit due to "H^?^ (Eq. [28|) : 

IPN precession of the outer orbit due to (Eq. [30|): 



and a new IPN interaction term between the two orbits, 
-Hj^™ (Eq. EJ), which introduces a new inclination and 
eccentricity dependent modulation (e.g., Fig. |3]). 

We compared the different timescales associated with 
the secular Newtonian and different IPN terms (see 
Fig. HI). If the timescales associated with the IPN effects 
are much shorter than the timescales associated with the 
eccentric Kozai-Lidov mechanism, i.e. the secular New- 
tonian timescales, the growth of the eccentricity in the 
inner orbit tends to be suppressed. We confirm that the 
excitation of the eccentricity is indeed suppressed for sys- 
tems where the Kozai-Lidov timcscale is many orders of 
magnitude longer than the IPN timescales. However, if 
the timescales of the IPN effects arc comparable to the 
secular Newtonian ones (see Fig. |4l) , we found two inter- 
esting regimes that present qualitatively different behav- 
ior. 

The first regime is where the IPN timescales are 
comparable but slight ly shorter than t he Newtonian 
Kozai-Lidov timcscale. iFord et all (|2000bf ). studying the 
PSR B1620— 26 triple system, noted that the inner ec- 
centricity may be greatly increased around some critical 
value of the outer SMA, due to the H^?? term and the 

octupole term. We extended this calculation by including 
all averaged IP N terms up to Oja '^) and the Newtonian 
octupole term (jNaoz et al.ll201l"b() . as well a s exploring 
a wide region of phase space. We confirmed iFord et al] 
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(|2000bD result and found a resonant-like behavior, where 
the inner orbital eccentricity is greatly increased com- 
pared to the Newtonian case. This behavior exists also 
when including all averaged IPN terms and for a wide 
range of mass ratios and orbital parameters. We param- 
eterized the location of the resonant peak in terms of the 
SMAs by defining a parameter, TZ in Eq. (|Tf)) . as the ratio 
of the leading-order IPN and secular Newtonian terms. 
This parameter depends on the ratio of the mass of the 
outer perturber to the total mass of the inner binary. 
The presence of the octupole term is important for the 
resonant IPN eccentricity excitation, which is most ap- 
parent in the examples with a small mutual inclination. 
For systems where either the inner or the outer binary 
shrinks, for example due to GW radiation-reaction, the 
triple may pass through this three-body IPN resonance. 
The amplitude and location of the resonance changes due 
to IPN terms as a function of TZ. We found that lower 
mutual inclinations in the prograde regime cause a wider 
peak (in terms of TZ), while a less massive outer body 
tends to produce wider and higher amplitude peaks. A 
detailed investigation of the properties of the resonance is 
beyond the scope of this paper, but could be the subject 
of future investigations. 

The second regime that exhibits qualitatively different 
behavior from that obtained with a quadrupole Newto- 
nian Kozai-Lidov treatment is when the quadrupolar sec- 
ular Newtonian timescales are shorter than the IPN ones 
and when the inner binary has comparable mass compo- 
nents. The eccentric Kozai-Lidov mechanism, neglecting 
IPN effects, is suppressed when mi — > 7712, since the 
outer orbit's potential is effectively quadrupolar. As we 
showed in this paper, IPN effects can break symmetry 
and excite eccentricity, triggering the eccentric Kozai- 
Lidov mechanism. As long as IPN precession occurs on 
a comparable timescale (or lower) than the Newtonian 
octupole precession, i.e, t^^^^^ £ ~ ^octJ the eccen- 
tric Kozai-Lidov mechanism will be triggered. 

Eccentricity excitations are particularly interesting 



in the context of pos si ble GW detections (iWenI 120031: 
Brown fc ZimmermanI l2010t lArmitage fc NataraianI 
2005H S esan al l2010() . If such excitations were not present, 
the frequency of the GWs emitted by the inner binary 
would be ty pically too low for detection with LIGO 
(see however lOTearv et all [20091: iKocsis fc Levh]l201lL 
for eccentric binaries which form in the LIGO band). 
However, if eccentricity is secularly excited through a 
three-body interaction, the frequency of the GWs is also 
increased during pericenter passage, thus bringing the 
signals into the detector's sensitivity band. Such large 
eccentricities would then lead to GW-driven inspiral and 
the eventual merger of binaries. Whether such eccentric 
signals can be detected or not will depend on how close 
such sources arc to Earth. But if detections are made 
with sufficiently high signal-to-noise ratio, then GWs 
could be used to measure the eccentricity of the inner 
binary, and thus, distinguish between different source 
populations. 
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APPENDIX 

TWO BODY SYSTEMS - EFFECTIVE ONE BODY 

GR pericenter precession has been studied in great detail and used to te st Einstein's theory in the Solar System, for 
example with observations of the perihelion precession of Mercury (e.g., iShapiro et al.l 119721 ). The simplest method 
to derive such precession is to consider test-particle motion in an effe ctive potentia.1, assu ming that GR introduces 
small corrections to N ewtonian dynamics and small eccentricities (e.g ., iMisner et al.|[l973l chapter 25 p. 668-6701^. 
The Hamiltonian (e.g. lArtemova et al.lll996l : IMiller fc HamiltonI 120021) is simply derived by integrating over the pre- 
cession rate. Althoug h this Hamiltonian leads to the correct ISGO location, if one uses the full expansion given in 
lArtemova et al.l (|1996l Eq. 4) , it is not equal to the IPN Hamiltonian. 

The purely orbital (non-spinning) 3PN Hamiltonian was derived in lJaranowski fc Schafeii (|1998l[200ll ) (in the center 
of mass frame, and after subtraction of the total rest- mass term). Here, we focus only on expansions to IPN order for 
a two body system (w ith masses mi and m2 and momenta pi and p2, respectively). The Hamiltonian is then (e.g., 
IBuonanno et al.ll2006[) : 

.^2body _ {ml + mDp^ P{3ml + 7mim2 -|- SniDp'^ k'^{p ■ r)^ ^ fc'*(mi + m2)^/^ 
8c^mi^m2 2c^Tnim2r 2c'^r^ 2(?r'^ ' 

where r is the radius vector between the two bodies, with magnitude r, the linear momentum of the effective one body 
problem is simply p = —pi = p2 and /i = TOim2/(mi -f- m2). Eliminating the short-period terms in the Hamiltonian, 



Note that the same pr ecession rate can be also derived directly 
from the IPN metric (e.g.. IMisner et al.|[T973 . chapter 40 p. 1100- 
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using the Von Zeipel transformation (for more details, see iBrouweHllQSOl ) for an orbit with SMA a and eccentricity e, 
the double average Hamiltonian is given by 



n 



which is the same as Eq. 



2body _ k^^^in (l5mi^ + 29mim2 + 15m2^) 3fc^TOim2(mi + 1712) 



IPN 



(Al) 



EQUATION OF MOTIONS FOR THE IPN INTERACTION 

Using the canonical relations [eqs. (|T5]) ]. we find the equations of motion for the interaction part of the IPN Hamil- 
tonian: 



dgi 
dt 



iPN(int) 4a^c2 (1 - 62^)^^^ (mi + 7712) 
aifc^mim2m3 



G2(8(toi + 7712) + 61113) COSitot 



(mi + m2) 



l_ei 
Gi 



(fmei (1-3 cos^ itot) + 18ei(mi + mim2 + ml) cos(2gi) sin^ Hot j 

aik^mim2m3 , 



GiG2(8(mi + m2) + Gms) sin it 



where 



+ 9ei(mi + mim2 + ml) cos(2(7i) sin(2itot) 

1 dfjYiei 



8(mi + m2) 
1 , cos itot 



3sin(2ltot)/mei 



G2 sin itot Gi sin itot 



fmei 



ei dt 



lOmi + 6mim2 — 10m2 



dg2 fc^ I G2 + Gi cos itot 

dt ~ 32a3c2(l - e2)3/2G2(ml + m2) \ Gi(mi+m2 



16GiG2(mi + m2)^4 



-|-6aifc^m2m3 cos itot (—/mei"^i — Ze\mi{m\ + mim2 + m\) cos(2(7i)) 



1 



mi + m2 

+ %e\{m\ + mim2 + m2) cos(2gi) sin^ itot) 
where /i4 = 4 (mi + TO2) + 3m3, 



32/14 cos itot (mi + m2) + 3aifc^mim2m3(/„iej - 3/mei cos^ itot 



dei 
~dt 



9aieiA/l — e\k'^mim2{m\ + mim2 + m\)m3 siv? itot sin(2(7i) 



iPN(int) 16a3c2(l - e2)3/2Li(mi + m2)2 

The change of the inner orbital angular momentum is simply 

%aie\k'^mim2{m\ + mim2 + m\)m3 sin^ itot sin(2(7i) 

IPN(int) 



dGi 



16a3c2(l-e2)3/2(mi +m2)2 



while for the outer orbit it is simply zero. Thus, 

dHi 



dt 



IPN(int) 



sin i2 dGi 
sin itot dt 



I.e., 



dH^ 



dt 



sini2 Qaie\k mim2{m\ + mim2 + m^ms sin itot sin(2(7i) 



iPN(int) sinitot 16a^c2(l - e|)3/2(mi + m2)2 

The inclinations evolve according to (cosii) = Hi/Gi — Gi/Gi cosii fe.g. lNaoz et al.|[2011bf ). and thus, 

9aieifc*miTO2(mi + mim2 + m'^)m3 sit? itot sin(2(7i) 

IPN(int) 



(cosii) 



16a3c2(l --e2)3/2(mi +m2)2 



1 / sin i2 



Gi V sinitot 



cos«i 



(Bl) 



(B2) 
(B3) 



(B4) 



(B5) 



(B6) 



(B7) 



(B8) 
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and since (00312) — H2/G2 — G2/G2Cosi2 and H2 ~ —Hi (e.g. iNaoz et al.|[2011bD we find 

•., sin 12 9aie^fc^mim2(mi + mim2 + m|)7n3 sin^ itot sin(25i) 

" G2 sin Hot 16a3c2(l-e2)3/2(mi+m2)2 ' ^ ^ 
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